import matplotlib.pyplot as plt
import numpy as np

# 三次贝塞尔曲面
def double3bezier(Pts):
    n = 100
    Mbe = np.array([[-1, 3,-3, 1],
                    [ 3,-6, 3, 0],
                    [-3, 3, 0, 0],
                    [ 1, 0, 0, 0]])
    P = np.transpose(Pts, [2,0,1]).reshape(1,1,3,4,4)
    arr = np.linspace(0, 1, n)
    arr = np.array([arr**3, arr**2, arr, np.ones((n,))])
    U = np.transpose(arr.reshape(4,1,n), [2,1,0]).reshape(1,n,1,1,4)
    W = np.transpose(arr.reshape(4,n,1), [1,0,2]).reshape(n,1,1,4,1)
    '''
        U  : (1, n, 1, 1, 4)
        W  : (n, 1, 1, 4, 1)
        P  : (1, 1, 3, 4, 4)
        Mbe: (1, 1, 1, 4, 4)
    '''
    points = U @ Mbe.reshape(1,1,1,4,4) @ P @ (Mbe.T).reshape(1,1,1,4,4) @ W
    return points.reshape(n,n,3)


if __name__ == '__main__':
    # 创建画布
    plt.rcParams["font.sans-serif"]=["SimHei"]  # 设置字体
    plt.rcParams["axes.unicode_minus"]=False    # “-”负号的显示
    fig = plt.figure()                          # 空的画布
    ax = fig.add_subplot(111, projection='3d')  # 创建3d绘图区域
    ax.tick_params(labelsize=8)                 # 刻度字体

    # 绘图
    Pts = np.array([[(0, 0, 0), (0, 4, 0), (0, 8,-3), (1, 12, 5)],
                    [(2, 0, 6), (2, 4, 0), (2, 8, 0), (2, 13, 7)],
                    [(4, 0, 0), (4, 4, 0), (4, 8, 3), (4, 15,10)],
                    [(6, 0, 0), (6, 4,-3), (6, 8, 0), (6, 11, 8)]])
    points = double3bezier(Pts)
    ax.plot_surface(points[:,:,0],points[:,:,1],points[:,:,2])
    ax.set_title('三阶贝塞尔曲面')
    for p in Pts:
        ax.plot(p[:,0],p[:,1],p[:,2])
    for p in np.transpose(Pts, [1,0,2]):
        ax.plot(p[:,0],p[:,1],p[:,2])

    # 参考线
    def axs(xx=(-1,1),yx=(-1,1),zx=(-1,1)): 
        color = ['r', 'g', 'b']
        for i,t in  enumerate((xx, yx, zx)):
            xs = [(0,0), (0,0), (0,0)]
            f = True if t[0]<0 and t[1]>0 else False
            if f:
                xs[i] = (t[0], 0)
                ax.plot(xs[0], xs[1], xs[2], linestyle=':', color=color[i])
                xs[i] = (0, t[1])
                ax.plot(xs[0], xs[1], xs[2], linestyle='-', color=color[i])
            else: 
                xs[i] = (t[0], t[1])
                ax.plot(xs[0], xs[1], xs[2], linestyle='-', color=color[i])
    axs((0,5),(0,5),zx=(0,5))
    ax.grid(b=None) # 隐藏网格

    # 显示窗口
    plt.show()

